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Prologue 


This TM consists of two chapters. Chapter I describes the development of a coupled, one-dimensional 
biogeochemical model using turbulence closure mixed layer dynamics. The model is applied to the 
Sargasso Sea at the BATS (Bermuda Atlantic Time Series) site and the results are compared with a previous 
model study in the same region described in NASA/TP-200 1-209991. 

Chapter II addresses the sensitivity of global sea-air C0 2 flux estimates to wind speed, temperature, and 
salinity. This is a revision of a previous analysis documented in NASA/TM -2002-211604. Sensitivity 
analyses of sea-air CO, flux to gas transfer algorithms, climatological wind speeds, sea surface temperature 
(SST) and salinity (SSS) were conducted for the global oceans and selected regional domains. Large 
uncertainties in the global sea-air flux estimates are identified due to different gas transfer algorithms, 
and seasonal wind speed, SST and SSS climatologies. The global annual sea-air flux ranges from -0.6 to 
-2.3 Gt/yr, depending on the combination of gas transfer algorithm and climatological wind speeds used. 
However, two of the gas transfer algorithms narrow the range of uncertainty for the estimated sea-air flux 
(-l.i a -1.4 Gt/yr). There are significant geographical differences in the distribution of oceanic sources 
and sinks of CO, originating from distinct regional patterns of the seasonal pC O, variability. Different 
combinations of SST and SSS global fields resulted in changes as large as 16% on the oceans global sea-air 
flux. This study shows that climatic changes in the global heat balance may significantly affect the global 
C0 2 flux via changes in SST. This study also highlights the potential impacts of changes in the oceanic 
hydrologic cycle on C0 2 fluxes, although salinity changes would be regional and could be either positive 
or negative. SSS changes in regions of high flux could be significant to the global CO, flux. 



Abstract 


This TM consists of two chapters. Chapter I describes the development of a coupled, one-dimensional 
biogeochemical model using turbulence closure mixed layer (TCMLM) dynamics. The model is applied 
to the Sargasso Sea at the BATS (Bermuda Atlantic Time Series) site and the results are compared 
with a previous model study in the same region described in NASA/TP-200 1-209991. The use of the 
TCMLM contributed to some improvements in the model simulation of chlorophyll, PAR, nitrate, 
phosphate, and oxygen, but most importantly, the current model achieved good agreement with the data 
with much more realistic background eddy diffusivity. However, off-line calculations of horizontal 
transport of biogeochemical properties revealed that one-dimensional dynamics can only provide a 
limited assessment of the nutrient and carbon balances at BATS. Future studies in the BATS region will 
require comprehensive three-dimensional field studies, combined with three-dimensional eddy resolving 
numerical experiments, to adequately quantify the impact of the local and remote forcing on ecosystem 
dynamics and carbon cycling. 

Chapter II addresses the sensitivity of global sea-air C0 2 flux estimates to wind speed, temperature, and 
salinity. Sensitivity analyses of sea-air C0 2 flux to wind speed climatologies, gas transfer algorithms, SSS 
and SST were conducted for the global oceans and regional domains. Large uncertainties in the global 
sea-air flux are identified, primarily due to the different gas transfer algorithms used. The sensitivity of 
the sea-air flux to SST and SSS is similar in magnitude to the effect of using different wind climatologies. 
Globally, the mean ocean uptake of C0 2 changes by 5 to 16%, depending upon the combination of SST 
and SSS used. 
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Chapter I 


Development of a Coupled, One Dimensional Biogeochemical Model 
Using Turbulence Closure Mixed Layer Dynamics and Application to 
the Sargasso Sea (BATS site) 


1.1. Background 

The mixed layer component of an existing one dimensional biogeochemical model 
[Signorini et al., 2001a; Signorini et al., 2001b], hereafter referred to as ECOID, was 
replaced with a turbulence closure mixed layer model [Mellor and Yamada, 1982]. 
ECOID has a bulk model of the mixed layer [Garwood, 1977] such that an explicit 
profile of the vertical eddy diffusivity (k v ), required by the biogeochemical conservation 
equations, cannot be provided by the model. Instead, an ad-hoc analytical fit scaled to 
surface wind stress, usually an exponential decay function that matches the mixed layer 
k v value with a much smaller background eddy diffusivity, was used to provide values of 
k v at each depth (Az=l m). The turbulence closure mixed layer model (TCMLM) has a 
vertical coordinate system that provides parameter values, including k v , at each time step 
and grid point. Also, the surface boundary layer is resolved more accurately by using a 
stretched vertical coordinate with higher resolution near the surface. Therefore, the 
coupling between the physical and biogeochemical processes can be achieved in a more 
robust and accurate manner. 


The new model, consisting of the TCMLM coupled to the biogeochemical conservation 
equations, was rewritten in Fortran 90/95 [Absoft, 1997] for a Linux machine. The 
following sections describe the details of the model development and its application to 
the BATS (Bermuda Atlantic Time Series) site. 

1.2. Mixed Layer Model 

The level 2-1/2 version of the Mellor- Yamada turbulent closure scheme was modified to 
account for solar flux divergence following the method of Simpson and Dickey [1981]. 
The method was further modified to account for spectral decomposition of the downward 
radiation. Details of the spectral method are described in the following section. The 
details of the closure scheme may be found in Mellor and Yamada [1982]; however, a 
brief description of the model will be included for completeness. Assuming hydrostatic 
approximation, the mean zonal and meridional momentum equations are 
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where u and v are the mean zonal and meridional velocity components, u g and v g are the 
corresponding geostrophic velocities, -u'w' and -vV are the vertical turbulent 
Reynold stresses, ju b is the background diffusivity, / is the Coriolis parameter and t is 

time. The last terms on the right side of (1) and (2) are the vertical advection of 
momentum. The vertical velocity w is pre-calculated for model forcing (see Section 5 for 
details). 

The heat conservation equation, which includes the divergence of the downward 
irradiance in spectral form, is given by 


dT 
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dm 

dz 
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where T is the mean temperature and c p is the specific heat of seawater. The first term is 

the local rate of change of heat, the second term is the vertical divergence of turbulent 
heat flux. The third term is the vertical advection of heat, and the fourth term is the pre- 
calculated horizontal heat divergence using climatological currents and temperature 
three-dimensional fields (see Section 5 for details). The right-hand side of (3) is the 
vertical divergence of the downward irradiance flux. The vertical distribution of 
downward irradiance is given by 

J -3000 eh f 3000 

1=280 7(/l,Z) = f=0 I= 2 80 7 0 ) ex p[a(^)^]^ dz (4) 

where 7(A,z)is the irradiance per unit wavelength for wavelength T (nm)at depth z, 
a(X) is the irradiance attenuation coefficient (m' 1 ) for wavelength X , and z is the vertical 
coordinate positive upward with origin at mean sea level. 

The salt conservation equation can be written as 

dS d TTTi dS dS — . 

— + — (w 5 + — ) + ^— + U cim • V5 clm - Q F (5) 

dt dz dz dz 


where S is the salinity and Q, represents the sources and sinks of freshwater. As in (3), 
the vertical advection and horizontal divergence of salt are also included. At the surface, 
Q, {Qi) = 5'(0)[£'-/ > ]/cfe(0); Q h =0at all other levels below the surface. Therefore, the 

salinity equation is coupled to atmospheric moisture divergence/convergence via 
evaporation minus precipitation ( E-P ). The water density, p(T,S,z), is calculated 
using an equation of state. 

Turbulent kinetic energy (TKE), 1/2 q 2 =1/2 («' 2 + v' 2 +w' 2 ) , and a turbulent length scale 
/ are computed. The level 2-1/2 version used here includes the diffusion of TKE, which 
appears to be important for low wind speed situations. The equation for TKE may be 
written as 
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( 6 ) 
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where the first term is the tendency of TKE, the second term is the diffusion of TKE, the 
first and second terms on the right side of (6) are the shear production, the third term is 
the buoyancy production and the last term ( D ) is the dissipation. 

The problem is closed by invoking hypotheses and assumptions regarding third moment 
terms [. Mellor , 1973]. Required empirical constants are determined from laboratory flow 
and are not adjusted. The details of the model are rather involved. However, it is 
possible to reduce the turbulent Reynolds stress and heat flux terms to the following 
expressions 
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Sm and Sh are stability functions; k m and k h are the eddy viscosity and eddy diffusion 

coefficients, respectively. The vertical diffusion coefficient for the biogeochemical 
model is obtained by adding the eddy diffusivity to the background diffusivity 
(k v =k h + p b ). The background eddy diffusivity (minimum value attained below the 

surface mixed layer) can be justified as a parameterization of processes not included in 
the model such as internal waves [ Chen and Annan , 2000]. Oceanographic estimates for 
the background diffusivity are in the range 10’ 5 - 10‘ 4 m 2 s’ 1 [ Christian et al., 1997; 
Mourn and Osborn, 1986; Sharpies et al., 2001; Siegel and Domaradzki, 1994]. A deep 
ocean estimate of an upper bound to the eddy coefficient for vertical diffusion is 10’ 4 m 2 s’ 1 
[Mourn and Osborn, 1986]. 

The flux Richardson number, which is simply the ratio of negative turbulent buoyancy 
production to shear production, is defined as 
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An important feature is the extinction of turbulence (Sm =Sh = 0) for R if > 0.21 [Mellor, 
1973]. 
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The numerical scheme uses a two time step implicit algorithm [Richtmyer and Morton, 
1967]. The vertical grid spacing is optimized by using a logarithmic distribution, which 
is ideal for the input of the penetrative component of solar irradiance. For the BATS 
application, there are 41 levels for a depth of 350 meters with the density of grid points 
being greatest near the surface. The time step is 30 minutes. 

The boundary conditions for the model (in addition to the surface value for Q h in (5)) 
include 

1 ) wind stress t ox , r oy : 


~ rjii 

T ox =pqlS M {R lf ) — 

2=0 

(11) 

T v =pql~S u (Rj)^ 

uZ 2= o 

(12) 


2) Turbulent heat flux H: 

(13) 

z=0 

The model is initialized with specified temperature, salinity, and horizontal velocity 
profiles. 

1.3. Irradiance Model 

The total (infrared plus visible) solar radiation is obtained using the Frouin model 
[Frouin et al., 1989], This model provides the total radiation and the photosynthetically 
available radiation (PAR). The infrared (I IR ) component is obtained by subtracting the 
PAR component from the total solar radiation. Using a spectral model for PAR [Gregg 
and Carder, 1990], the spectral PAR component I rAK (A) can be determined. The infrared 

component (for mixed layer model only) and the PAR component (mixed layer and 
biogeochemical model) of the penetrating irradiance are obtained from 

I, R (z) = I m (z-Az)exp[-a m Az] (14) 

I PAR (*,*) = I (A, z - Az) exp[-(a w (A) + a ph ( A))Az] (15) 

where a IR ( 3.75 m' l )is the attenuation coefficient for infrared radiation, and a w (A) and 
a ph (X) are the wavelength-dependent light attenuation coefficients for water and 


H = pc p qlS H (R lf ) — 
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phytoplankton, respectively. The model has two options for the water attenuation 
coefficients based on Baker and Smith [1982] and Pope and Fry [1997]. Two options for 
chlorophyll-dependent attenuation coefficients are available in the model, those from 
Morel [1988] and those from Dupouy et al. [1997]. Figure 1-1 shows spectra for water 
absorption and chlorophyll absorption coefficients using the available model options. 
Another option is the inclusion of the dissolved matter attenuation coefficients, a im (A) . 

The dissolved matter attenuation coefficients are calculated by applying the IOP (inherent 
optical properties) model of Garver and Siegel [1997], which uses water leaving 
radiances from 6 SeaWiFS bands as input (level 3 binned monthly composites). The IOP 
model calculates the attenuation coefficient due to dissolved matter for the 443nm 
wavelength. The attenuation coefficients for other wavelengths are obtained from 

W = « dm (443) exp [S(T - 443)] ( 1 6) 

where the exponential decay constant, S , is chosen to be 0.02061 . 

A correction is applied to the downward irradiance pathway to account for seawater light 
refraction following Snell’s law. After some algebraic manipulations, the correction is 
applied to Azas follows 


Az' = 


Az 

cos/? 


(17) 


P = a sin 



(18) 


where a. is the solar zenith angle, and n s is the seawater refraction coefficient which is 

expressed as a function of salinity and temperature adapted from Table 3.12 of Neumann 
and Pierson [1966] as 


n s = 1 (T 6 (285.77 - 1 5.65 T + 1 97.67 S) + 1 .333338 ( 1 9) 


1.4. Biogeochemical Model 

The biogeochemical model is an upgraded version of the Signorini et al. [2001a] model, 
which includes additional conservation equations for diazotrophs and detritus. The 
addition of diazotrophs and detritus in the model follow the method of Fennel et al. 
[2002]. Figure 1-2 shows a diagram illustrating the model components and their 
couplings. The governing equations for the biogeochemical model are 
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Figure 1-2. Flowchart of model components showing the interaction between the nitrogen, 
phosphorus, and carbon pools. 
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The subscripts phy, zoo, nf, det, and elm refer to phytoplankton, zooplankton, nitrogen 
fixers, detritus, and climatology, respectively. Table 1-1 defines the model state variables 
and Table 1-2 provides the definition of the model parameters and values used. 


Table 1-1. Model state variables. 



Unit 

Variable 


no 3 


Dissolved inorganic nitrogen 


P0 4 

id kit 

Dissolved inorganic phosphorus 


ESH 


Ammonium 



pmol kg" 1 

Dissolved oxygen 



WBSSESKi 

Dissolved inorganic carbon 


NF 


Nitrogen-fixing phytoplankton 


P 


Other phytoplankton 


Z 


Zooplankton 


DON 


Dissolved organic nitrogen 


DOP 


Dissolved organic phosphorus 


DOC 

nmol kg" 1 

Dissolved organic carbon 


det N 


Detrital particulate nitrogen 


detp 

mmol P m° 

Detrital particulate phosphate 



1.5. Application to the BATS site 

The model was applied to the Bermuda Atlantic Time Series (BATS) site and the results 

compared with a previous version of the biogeochemical model [Signor ini et al . , 2001a] 

with off-line physical forcing (temperature, salinity, and mixed layer depth) derived from 

CTD casts. The major additions and modifications to the current model are 

• The physical forcing parameters, temperature, salinity, mixed layer depth, and 
vertical eddy diffusivity, are now computed using the TCMLM 

• The Ekman pumping velocity was derived from daily NCEP winds. The previous 
model used Ekman pumping velocity derived from monthly winds. 

• The horizontal divergence of heat, salt, nitrate, phosphate are estimated from seasonal 
climatology of currents (geostrophic plus Ekman drift), nitrate, and phosphate (see 
Figure 1-4) 

• The horizontal divergence of dissolved inorganic carbon (DIC) is estimated from 
seasonal climatology of currents and the DIC algorithm of Lee et al. , which is 
dependent on temperature, salinity, apparent oxygen utilization (AOU), and silicate, 
also obtained from seasonal climatology (see Figure 1-4). 

• A dilution term (salinity-dependent) was introduced for the DIC conservation 
equation (3 1 ) 
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Table 1-2. Model Parameters. 

Symbol I Value I Unit 



mmol N: mmol P 




Dimensionless 


Dimensionless 


mmol N d 


WnT 


Parameter 

Phytoplankton respiration rate 


Phytoplankton mortality rate 
Half-saturation concentration for N0 3 uptake 


Half-saturation concentration for NH 4 uptake 


Half-saturation concentration for P0 4 uptake 


Stoichiometry of phytoplankton and zooplankton 


Maximum growth rate of phytoplankton 


Initial P-I slope of phytoplankton 


Diazotrophic respiration rate 


Diazotrophic mortality rate 


Maximum growth rate of diazotrophs 


Half-saturation concentration for P0 4 uptake by diazotrophs 


Initial P-I slope of diazotrophs 


Critical temperature for diazotroph growth 


Critical wind stress for diazotroph growth 


Stoichiometry of diazotrophs 


Maximum grazing rate of zooplankton 


Half-saturation concentration for ingestion by zooplankton 


Zooplankton assimilation efficiency 


Zooplankton excretion rate 


Zooplankton mortality rate 


Remineralization rate of detritus 


Sinking velocity of detritus 


Remineralization rate of DON 


Remineralization rate of DOP 


Remineralization rate of DOC 


Fraction of dead phytopalnkton converted to ammonium 


Fraction of dead zoopalnkton converted to ammonium 


Fraction of dead phytopalnkton converted to DOM 


Fraction of dead zoopalnkton converted to DOM 


Maximum rate of ammonium nitrification 


Maximum light inhibition dosage for nitrification 


Half-saturation light dosage for nitrification photoinhibition 




















































• A correction for light penetration refraction was added based on Snell’s law. 

• The chlorophyll-dependent a ph (A) from Morel [1988] was added as an option. This 
option was used in the current simulation. 

• An option for light absorption by dissolved matter derived from an IOP model and 
SeaWiFS water leaving radiances was developed (not used in current simulation). 

• Values for DIC and O 2 at the bottom boundary are no longer constant; they were 
derived off-line from the BATS data time series. The NO 3 and PO 4 bottom boundary 
values are also obtained similarly, as in the previous model simulation. 

• New equations were added to the biogeochemical model for determining the 
concentrations of diazotrophs, and detritus (nitrogen and phosphorus components). 


The following subsections describe the model configuration, forcing, validation, and 
comparison of results with the previous model. 

Data Sources and Model Forcing 

The ecosystem model is physically forced by temperature, salinity, and mixed layer depth 
from the TCMLM. The BATS biogeochemical data were used for model forcing and 
validation. The data are discussed in several publications [Bates et al., 1996; Michaels 
and Knap, 1996; Michaels et al., 1994]. The observed PAR profiles for validation were 
obtained from the Bermuda Bio-Optics Program (BBOP) database. The 1994-1996 
seawater and atmospheric pC 02 data originate from the carbon dioxide measurements 
obtained from the R/V Weatherbird II in the Sargasso Sea [Bates et al., 1998]. Most of 
these data sets were obtained interactively via the web-based data extractor maintained by 
BBSR (Bermuda Biological Station for Research, Inc.; 
www.bbsr.edu/cintoo/bats/bats.html). The eddy-induced vertical velocity time series was 
calculated (see Section 3) from the merged TOPEX/Poseidon-European Remote Sensing 
Satellite (ERS), 1/3°, 7-day sea-surface height anomaly (SSHA) grids [AVISO, 1997]. 
The SSHA data were obtained via anonymous ftp from the AVISO data distribution 
server (ftp://ftp.cls.fr/pub/oceano/AVISO). Three-hourly wind speed data for the gas 
transfer formulation, and for mixed layer model forcing, were obtained from the Bermuda 
Airport weather station. The gridded winds to calculate the Ekman pumping vertical 
velocity, and cloud cover data to calculate surface PAR, and meteorological parameters 
(precipitation, air temperature, and humidity) to calculate surface heat and freshwater 
fluxes for the TCMLM, were obtained from the NCEP/NCAR Reanalysis Project 
[Kalnay et al., 1996] daily products. 

The BATS region is characterized by energetic mesoscale eddies that provide vertical 
transport of nutrients to the euphotic zone via eddy pumping [McGillicuddy et al., 1999; 
McGillicuddy et al., 1998a; Siegel et al., 1999]. The vertical velocity variability in the 
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BATS study area is affected primarily by vertical displacements of the thermocline 
imparted by the migration of mesoscale eddies through the area [McGillicuddy et al., 
1998b; Siegel et al., 1999]. This eddy pumping mechanism has been demonstrated to 
have a significant impact on the nutrient budget at BATS [Siegel et al., 1999]. However, 
Oschlies [2002], using a three-dimensional eddy-resolving (1/9° resolution) model, 
concluded that eddy pumping has a much more limited role in nitrate supply at BATS 
than reported by previous studies. We obtained the eddy pumping vertical velocity 
( w eddy ) time series for the BATS location by taking the first derivative of the SSHA and 

scaling it to the pycnocline displacement using the scale factor of 4 meters of thermocline 
displacement for each cm of SSHA. The scale factor originates from Siegel et al. [1999]. 

Figure 1-3 shows time series of geostrophic currents, SSHA, and vertical velocity for the 
available time period of 21 October 1992 through 13 February 2002. An uninterrupted 
time series of eddy vertical velocity, covering the period of 1990-1998, was constructed 
by filling the missing points (1 January 1990 - 20 October 1992) with data from adjacent 
years. The model was run for the period 1 January 1990 - 31 December 1998, but the 
first two years of simulation were discarded as they served as spin-up years. The vertical 
velocity profile, w(z), is parameterized as, 


w (z) = Weddy ^, for 2 < 200m 
w(z) = w ecjdy , for 200w < z < 350w 


(48) 


The 200 m transition depth follows the [Siegel et al., 1999] empirical parameterization. 
Below 200 meters isopycnal displacements are related to SSHA as a 4-m isopleth 
displacement for each 1-cm of sea level height anomaly, while above 200 meters the 
displacement is modeled as a decreasing linear profile that extrapolates to zero at the sea 
surface. The weekly eddy-induced vertical velocity ranges from -5.3 to +5.9 m d’ 1 (solid 
line), while the daily Ekman pumping ranges from -0.95 to +0.62 m d' 1 (dotted line). The 
two components of the vertical velocity were added to force the model. The kinematic 
approach for eddy pumping derived by Siegel et al. [1999] requires the determination of a 
minimum SLA for a given uplifting event and assumes that all the nitrate brought to the 
euphotic zone is instantaneously used by phytoplankton. To meet these requirements, the 
eddy pumping events to force the model are obtained by specifying thresholds for SLA 
and w. The eddy pumping events are “turned on” whenever the SLA is negative 
(pycnocline uplifting) and the vertical velocity is positive. As an option, the thresholds 
for SLA and w can be turned off in the model allowing for upwelling (positive w) and 
downwelling (negative w) events to occur with no constraints. A comparison of the two 
methods will be offered in the forthcoming section discussing sensitivity tests. 

The geostrophic velocity components used in the TCMLM momentum equations are pre- 
calculated using the thermal wind equations at grid points surrounding the BATS site. 
The vertical profiles of geostrophic velocity are obtained from a combination of 
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TP/E RS- Derived Geostrophic Current, MSLA, and w at BATS 
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Figure 1-3. Time series of geostrophic currents (top), mean sea level anomaly and eddy pumping 
velocity (middle), and Ekman vertcal velocity (botom). The white line in the bottom tier is the 
low-pass-filtered Ekman velocity. 
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hydrographic data [ Conkright et al., 1998] from seasonal climatologies (baroclinic 
component) and SSHA (barotropic eddy component). The horizontal divergence of T, S, 
NO3, PO4, and TCO 2 (Figure 1-4) were calculated off-line using seasonal climatological 
data. When applied to the model, the horizontal divergence of PO4 derived from the 
seasonal climatology caused large discrepancies between the model and observed PO4 
profiles. The cause of these discrepancies is not clear. An alternate solution was adopted 
by obtaining a phosphate versus nitrate linear regression using the BATS data. Nitrate 
and phosphate are strongly correlated at BATS (r 2 =0.96) given by 

P 0 4 = 0.025449 + 0.052924NO 3 (49) 

The seasonal profiles of TCO 2 were obtained using the method of [Lee et al . , 1997]. 
Specifically, 


TC0 2 =1818 + 8.95S-2.33r-0.514AOU+0.54SiO 2 (50) 

where T, S, AOU and SiC >2 are the World Ocean Atlas 1994 seasonal temperature (7), 
salinity (5), apparent oxygen utilization (AOU), and silicate (Si 02 ). Data from the World 
Ocean Atlas 1994 were used instead of the World Ocean Atlas 1998 because the latter 
data set does not include AOU. To calculate the net horizontal divergence, the sum of 
(udC/dx + vdC/dy), where C represents the variables T, S, NO 3 , PO4, was computed on a 
2x2-degree grid surrounding the BATS site. The seasonal horizontal velocity 
components ( it, v, profiles) were obtained from the sum of geostrophic component (from 
hydrographic data) and Ekman drift derived from climatological wind stress components 
[< Oberhuber , 1988] and Ekman dynamics. 

The seasonal profiles of PAR, eddy diffusivity, vertical velocity, and u (east) and v 
(north) components of the velocity are shown in Figure 1-5. Except for u and v, the 
parameters shown in Figure 1-5 are used to force the biogeochemical model. 

Sensitivity Tests 

The model was run for the period 1990-1998, with the first two years of simulation 
discarded as spin-up years. The period of 1992-1998 was chosen because it corresponds 
to the period when most of the observed biogeochemical parameters are concurrently 
available at BATS. Numerous sensitivity tests were performed with the mixed layer and 
biogeochemical models to choose the best set of forcing approaches and parameter values 
that optimize model performance. Table 1-3 summarizes the major results. One of the 
key findings was the model’s inability to provide deep enough convective mixing during 
winter-spring to provide the required supply of nutrients to the euphotic zone. This 
model deficiency was concluded from five different runs that highlight the effect of 
vertical mixing on nutrient supply and primary production. The five runs are configured 
as follows. In run 1 the eddy diffusivity profile generated by the TCMLM is interpolated 
at the 1 -meter grid resolution and applied to the ecosystem model. The nitrate and 
phosphate profiles are computed prognostically without any constraints. The 
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Table 1-3. Summary of major model parameter sensitivity tests and outcomes. 


Parameters 

Sensitivity 

Comments 

Number of layers for the 
stretched vertical 
coordinate of the mixed 
layer model. 

The winter mixed layer depth, and 
therefore SST, is very sensitive to 
the number of layers used. 

This sensitivity remains illusive. 
The best SST values were obtained 
with 40 layers, which were 
hardwired in the original MLM 
code. The use of uniform layer 
spacing generates instabilities. The 
model requires a log distribution of 
grid points near the surface to 
resolve the air-sea boundary. 

Pre-calculated horizontal 
divergence of nitrate and 
phosphate from seasonal 
climatology. 

The resulting nitrate and phosphate 
profiles showed unrealistic 
behavior. Large nutrient influx in 
summer causes the spring bloom to 
shift to summer-fall. 

The horizontal divergence terms 
were turned off. 

Eddy pumping. 

Primary production is very 
sensitive to eddy pumping. Two 
approaches were used for w, eddy 
upwelling only {Siegel et al . , 1999) 
and unconstrained w (+/- values 
admitted). The two approaches 
bracket the ideal eddy pumping 
effect. 

A revised approach was 
implemented following Siegel et al 
( 1 999). The eddy pumping 
velocity is applied for uplifting 
events only (SLA<0 and w >0). 

Vertical eddy diffusivity 

The shape of the nutrient profiles 
and primary production are 
sensitive to this mixing parameter. 

Typical values for vertical eddy 
diffusivity are of order of 10' 5 m 2 /s 
{Gregg, 1987; Ledwell, 1993). An 
upper bound of 10' 4 m 2 s' 1 is given 
by Munk ( 1 966). A value of 5x 1 O’ 5 
m 2 s" 5 was chosen. 

Sea surface salinity. 

Most influential on the surface 
alkalinity and, therefore, an 
accurate estimate of surface TC0 2 . 

Surface salinity was strongly 
relaxed (order of 1 hour relaxation 
time scale) to the observed values. 

TC0 2 dilution. 

A function of salinity changes. 

The mean dilution effect on carbon 
flux is about 1/3 of the sea-air flux. 

Half saturation coefficients 
for nitrate and phosphate 
uptake. 

The magnitude and ratio of these 
two parameters determines the 
phytoplankton growth and the 
stability of nitrate and phosphate 
concentrations in the euphotic zone. 

The value for the half saturation 
coefficient for nitrate uptake is the 
same used in the previous model. 
The corresponding value for 
phosphate was reduced to 1/10 of 
its original value to enable a stable 
nitrate concentration and 
reasonable growth. 

Nitrogen Fixation 

The diazotroph growth is modeled 
as a function of phosphate 
concentration, light, sea-surface 
roughness (i.e., wind stress), and 
temperature {Fennel et al . , 2002). 

Parameter values for diazotroph 
equation 21 are shown in Table 1.2. 
The predicted 1992-1998 annual 
average input to new production by 
N2 fixation is 0.0035 mol N/m 2 /yr. 
This value is close to the 0.004 mol 
N/m 2 /yr reported by Orcutt et al 
[20011 at BATS. ‘ I 
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Figure 1-4. Seasonal profiles of nitrate, phosphate, dissolved inorganic carbon, temperature, and 
salinity horizontal divergences derived from seasonal climatologies of currents, nutrients, carbon, 
and hydrography. 
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Figure 1-5. Seasonal profiles of PAR, vertical eddy diffusion, vertical velocity, and u (solid) and v 
(dashed) components of horizontal velocity. The black dots superimposed on the PAR profiles are 
the data points, with the dashed line representing the mean profiles obtained from seasonal 
averages of the data. 
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eddy-induced vertical velocity is also free of constraints (upwelling/downwelling mode). 
Run 2 is similar to run 1, except that the eddy-upwelling mode is enabled (upwelling 
only, following Siegel et al. [1999]) and maintained for the other two subsequent runs. 
Run 3 is similar to run 2, except that the nitrate and phosphate concentrations below 
100m are relaxed to their mean climatological values with a relaxation time scale of 1 
month. Run 4 is configured with the nutrient relaxation to climatology and an ad-hoc 
reconfiguration of the eddy diffusivity profile. The eddy diffusivity profile is 
reconfigured as follows. First, the mixed layer depth (D m id) was calculated by finding the 
first depth where ae(D m id) - cre(O) = aj AT; aj is the coefficient of thermal expansion at 
sea surface conditions and AT is chosen to be 0.5°C [Siegel et al., 1995]. The potential 
density (ae) is obtained from the temperature and salinity profiles generated by the 
TCMLM. This criterion yields mixed layer depths in winter-spring that are ~50m deeper 
than those provided by the criterion used in the TCMLM code, which is based on finding 
the first depth where T(D m id) - T(0) = 0.2°C. The original TCMLM eddy diffusivity 
profile has a local maximum within the mixed layer, being zero at the surface and the 
background value below the mixed layer (Figure 1-5). The reconfigured eddy diffusivity 
profile assumes a constant value within the mixed layer equal to the maximum value 
provided by the TCMLM, similar to the profile assumed in McClain et al. [1996], The 
background diffusivity is not changed from its original value. A 50-point (50m) running 
boxcar filter is applied to the profiles for smoothing the transition between the mixed 
layer and the background values. Finally, run 5 was forced with TCMLM eddy 
diffusivity, upwelling/downwelling mode, and nutrient relaxation to climatology. 

Figures 1-6 through 1-9 show the seasonal profiles of primary production, chlorophyll, 
nitrate, and phosphate concentrations obtained from runs 1, 2, 3, and 4 respectively. The 
seasonal profiles from run 5 are not shown because the profiles are very similar to those 
from run 3. However, more quantitative results of all five runs will be discussed later in 
this section. The solid lines are the model profiles while the dashed lines are the seasonal 
means derived from the data (black dots). Figure 1-6 shows the profiles from run 1 . The 
agreement with the data is poor because the net vertical nutrient fluxes into the euphotic 
zone are too low to maintain the biological demand. It is concluded that the altimeter- 
derived vertical velocity time series provides too much downwelling and therefore does 
not promote the required vertical nutrient flux. Figure 1-7 shows the profiles from run 2, 
which was configured with the eddy-upwelling only scheme turned on. Figure 1-7 shows 
an overall agreement with the data, with the exception of the summer-fall primary 
production which is underestimated by the model, and a much sharper nutricline when 
compared to the observed nitrate and phosphate profiles. The summer-fall primary 
production is under-predicted by the model because the surface phosphate concentration 
is entirely exhausted during those seasons. One possible reason for this discrepancy is 
the omission of horizontal advection of nutrients in the model. Figure 1-4 shows that both 
nitrate and phosphate divergences for summer and fall have positive values (influx) in the 
surface layer and negative values (outflux) below the surface. This vertical structure 
supports surface enrichment and deep depletion of nutrients, which could mitigate the 
surface phosphate deficiency and the sharper nutricline predicted by the model. 
However, the magnitude of the nutrient divergence calculated from climatological fields 
may be too large since our tests showed that they cause excessive accumulation near the 
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Figure 1-6. Seasonal profiles of primary production, chlorophyll, nitrate, and phosphate from run 
(1). The solid lines are model simulations while the black dots and dashed lines are observations. 
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Figure 1-7. Seasonal profiles of primary production, chlorophyll, nitrate, and phosphate from run (2). 
The solid lines are model simulations while the black dots and dashed lines are observations. 
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Figure 1-8. Seasonal profiles of primary production, chlorophyll, nitrate, and phosphate from run (3). 
The solid lines are model simulations while the black dots and dashed lines are the observations. 
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Figure 1-9. Seasonal profiles of primary production, chlorophyll, nitrate, and phosphate from run (4). 
The solid lines are model simulations while the black dots and dashed lines are observations. 
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surface and too much depletion at depth. In addition, a comparison between runs 1 and 2 
shows the importance of the eddy pumping effect on the vertical flux of nutrients. These 
runs indicate that two extremes in the magnitude of the vertical flux of nutrients are 
achieved by the use of the two approaches for vertical eddy velocity. The adaptation of 
the Siegel et al. [1999] method for eddy pumping into the model provides reasonable 
results (although the nutricline is unrealistically sharp), while the use of unconstrained 
(positive and negative values) eddy vertical velocity does not provide sufficient pumping 
of nutrients into the euphotic zone. Short of developing an ad hoc vertical velocity 
approach tuned for the BATS site, the appropriate parameterization of eddy pumping in 
the model remains an elusive task. However, note that the vertical diffusivity is weakest 
in summer and fall and, despite the strong subsurface gradients of nitrate and phosphate 
being too high, the nutrient flux to the surface remains too low. This seems to indicate an 
inadequate penetration of mixing, demonstrated by the results of run 4, which will be 
analyzed in the forthcoming discussion. 

Figure 1-8 shows the profiles from run 3, which includes a relaxation of the nitrate and 
phosphate concentrations to the climatological profiles below 1 00m. This relaxation is a 
proxy to the nutrient horizontal divergence that could not be accurately derived from the 
climatological data. As shown in Figure 1-8, the nutrient profiles simulated in this run are 
in better agreement with the data. However, the model primary production is low in all 
seasons. This result leads to the conclusion that mixing in the model is too shallow 
during winter-spring and therefore limits nutrient supply. The reconfiguration of the 
eddy diffusivity profile for run 4 was designed to test the hypothesis of shallow mixing. 
Figure 1-9 shows the seasonal profiles obtained from run 4. All model profiles agree 
well with the data, except for the fall primary production, which is underestimated by the 
model but not as much as in the two previous runs. Although the reconfiguration of the 
eddy diffusivity provides the best parameterization for vertical mixing, it violates 
conservation of TKE in the TCMLM. As shown in equation (6), changing the eddy 
diffusion profile imparts in a change of buoyancy production (BP), as shown by 

BP = -S=w V = -=*,— (51) 

p p dz 

As indicated in equation (51), a change in eddy diffusivity {h) results in a change of BP 
and, therefore, a violation of TKE conservation. In view of this, even though the results 
of run 4 show a significant improvement on the performance of the ecosystem model, and 
they highlight a deficiency of the TCMLM, the procedure for reconfiguring the eddy 
diffusivity profile it is not internally consistent. Therefore, run 4 was not used in the final 
analysis. Instead, the results from run 2 will be used in the forthcoming analyses because 
no constraints were imposed on the nutrient computation and the eddy diffusivity 
obtained from the TCMLM is internally consistent. 

Nitrate Balance and the Relative Importance of Vertical Mixing and Advection 

The model nitrate balance in the upper 100m for the period of 1992-1997, during which 
there are available in situ concurrent observations of nutrients, primary production, and 
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chlorophyll concentration, is summarized in Table 1-4. The nitrate balance is calculated 
for the five different sensitivity runs, based on the model’s nitrate conservation equation 
(equation 25). Values from the literature are tabulated for available parameters. These 
include results from the eddy-resolving model of Oschlies [2002], and the statistical 
method of Siegel et al. [1999] based on altimeter data. 

The goal of this analysis is to quantify the relative roles of vertical mixing and vertical 
advection from mesoscale eddies on the nutrient supply to the euphotic zone. Table I-3b 
tabulates the nitrate flux for each term of equation 25 (nitrate uptake, ammonium 
remineralization, nitrogen-fixers (NF) respiration, vertical advection, and vertical mixing) 
for each of the five runs. The last column shows the ratio of mixing to advection. The 
NF respiration is relatively negligible. In order of importance, the main supply of nitrate 
is partitioned among vertical mixing, vertical advection, ammonium remineralization, and 
NF respiration. The ratio of mixing to advection varies widely from run to run. For the 
two runs that match the reported estimate of nitrate uptake of 0.5 1 5 mol N m' 2 yr'' [Siegel 
et al., 1999] for the Sargasso Sea more closely, runs 2 (0.592 mol N m' 2 yr'‘ ) and 4 (0.432 
mol N m‘ 2 yr'' ), the mixing to advection ratios are 1.2 and 4.5, respectively. Run 2 was 
forced with vertical velocity obtained with the upwelling-only mode and no relaxation to 
climatological nitrate concentration, while run 4 was conducted with relaxation to the 
climatological nitrate profile below 100m and the use of the ad-hoc eddy diffusivity 
profile, which provides more robust vertical mixing. The equivalent ratio derived from 
values of mixing and advection reported by Siegel et al. [1999] is 0.9, a value very close 
to the ratio obtained with run 2. The mixing to advection ratio derived from the eddy- 
resolving model of Oschlies [2002] is 4.4, a value very close to the ratio obtained from 
run 4 (4.5). However, the model of Oschlies [2002], despite its realistic rendition of the 
mesoscale eddy field, predicts a nitrate uptake (0. 1 6 mol N m' 2 yr- 1 ) for the BATS region 
that is in the low end of available estimates. Oschlies [2002] applied the statistical 
method of Siegel et al. [1999] to his model output at BATS and obtained a value for the 
vertical advection from mesoscale eddies which is much higher (0. 1 80 mol N m' 2 yr'' ) 
than the value predicted directly from his model (0.027 mol N m‘ 2 yr'' ). These available 
estimates reveal that there is a large uncertainty in the role of mixing and mesoscale 
eddies. One of the main reasons for this uncertainty stems from our inability to 
appropriately parameterize mixing and advection in our models. Despite the good 
agreement between model and observed SST and MLD, it was not possible to obtain a 
combination of realistic nutrient profiles and adequate supply of nutrients to the euphotic 
zone. A key element to resolve this problem is the design of field experiments at the time 
series sites (i.e., BATS and HOT) that will allow the determination of upper ocean (0- 
300m) eddy diffusivity profiles. More specifically, the shape of the eddy diffusivity 
profile within the mixed layer itself, and within the layer separating the base of the mixed 
layer from the more quiescent ocean below, must be determined to appropriately model 
vertical mixing. This is especially important for models with high vertical resolution 
such as the one used in this TM. The uncertainty in determining the shape and magnitude 
of the vertical velocity profiles derived from satellite altimetry is also a problem. These 
uncertainties will remain unresolved until new model parameterizations and field 
experiments can improve the realism of the physical forcing. 
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Biomass, Phytoplankton Nitrogen Production, and Variability of Phytoplankton Growth 

A common result of the sensitivity runs is that for most runs (except for run 1) the 
predicted biomass concentration matched the observations despite very low predicted 
primary production in 3 out of the five runs. To further investigate this apparent 
inconsistency, an analysis of the upper- 100m phytoplankton nitrogen flux (1992-1997 
mean) was conducted based on the balance of the source and sink terms in equation 20. 
The results are shown in Table 1-5. The phytoplankton nitrogen production is balanced, 
in order of importance, by zooplankton grazing, phytoplankton respiration and mortality, 
and physical transport by vertical advection and mixing. The physical transport is 
relatively small in all runs. The biomass and primary production are tabulated for all 
runs. The derived phytoplankton nitrogen production and phytoplankton biomass, are 
also shown in Table 1-5 (see footnotes in Table 1-5 for explanation on how phytoplankton 
nitrogen production and biomass were derived). 

The main result from this analysis is that there is a large variability of phytoplankton 
nitrogen production and grazing while the mortality, respiration, and biomass do not 
significantly change relative to growth and grazing. This indicates that the model 
biomass is under strong grazing control. This result is consistent with the behavior of the 
grazing, mortality, and respiration terms in equation 20. The zooplankton grazing term 
(equation 35) is proportional to the square of the biomass, while respiration and mortality 
(equation 20) are a linear function of the biomass. Therefore, it is possible to have low 
growth and reasonable values of biomass in the model, such that the model results must 
be evaluated for both biomass and primary production to accomplish a stringent 
validation test. For instance, run 2 is the only run that provides a combined reasonable 
match of new production (0.59 mol N m' 2 yr' 1 ), total phytoplankton nitrogen production 
(1.84 mol N nr 2 yr'' ), and biomass (14.8 mol N m' 2 ) with the available observed values. 
The corresponding observed values for the BATS region are 0.52 molNm‘ 2 yr ', 1.87 
mol N m‘ 2 yr'' , and 13.6 mol N m' 2 , respectively. 

An additional test was conducted to verify that the model formulations for the sources 
and sinks of phytoplankton nitrogen are appropriate for the BATS region. This was 
achieved by calculating the ratio of the standard deviation of the surface phytoplankton 
growth (/a in d 1 ) to the standard deviation of surface chl-a concentration (mg m‘ 3 ) for the 
model and observations. The standard deviations of both model and observations are 
based on the period 1992-1997. We calculate the phytoplankton growth from observed 
primary production (PP) and chl-a concentration at BATS using the formula 
//(d 1 ) = (PP chl:N)/(12 C:N chl-a) . The fixed C:N ratio (6.625) and the surface chl:N 
ratio (0.8) are the same values used in the model. The constant 12 is the molecular 
weight of carbon. The ratio of surface growth to surface chl-a standard deviations using 
values from run 2, which provided the best combined match of new production, 
phytoplankton nitrogen production, and biomass with the equivalent parameters derived 
from the BATS data, is 4.1 d' 1 (mg chl-a nr 1 )' 1 . The equivalent ratio derived from the 
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BATS data is 4.5 d' 1 (mg chl-a m' 3 )' 1 , which is 9.8% higher than for run 2. The similarity 
between model and data ratios shows that the balance between phytoplankton growth and 
biomass is accurately represented by the model dynamics. 


Table 1-4. Model nitrate balance in the upper 100m (1992-1997) from five different runs. Runs 1,2,3, and 
5 use k v from the TCMLM. Run 4 uses k v from the ad-hoc approach described in the text. Runs 1 and 5 are 
forced with unconstrained vertical velocity (upwelling/downwelling mode), while all the other runs are 
forced with upwelling-only mode {Siegel et al ., 1999). Relaxation to nitrate and phosphate climatology 
below 100m was imposed only for runs 3,4, and 5. All units are in mol N m" 2 yr' 1 . Values from the 
literature are tabulated for available parameters. 


Source 

no 3 

nh 4 

NF 

Vertical 

Vertical 

Mixing 


Uptake 

Rem. 

Respiration 

Advection 

Mixing 

Advection 

Run 1 


0.020 

0.00040 

0.006 

0.031 

5.2 

Run 2 

-0.592 

0.029 

0.00006 

0.251 


1.2 

Run 3 

-0.224 

0.026 




1.5 

Run 4 

-0.432 

0.030 



ixesihb 

4.5 

RunS 

-0.171 

0.020 

0.00003 

0.053 

0.098 

1.8 

Oschlies , 2002 

-0.16* 



0.027 

0.120 

4.4 

Oschlies , 2002 




0.180*** 



Siegel et al ., 1 999 

-0.515** 



0.240 

0.215 

0.9 


* Balance of 0.013 mol N m" 2 yr 1 due to other processes not reported in paper. 

** Includes Ekman horizontal flux and atmospheric deposition (0.06 mol N m" 2 yr 1 ). 

*** Application of statistical method of Siegel et al. [1999] to 3D model output by Oschlies [2002]. 


Table 1-5. Model phytoplankton nitrogen balance in the upper 100m (1992-1997) resulting from the five 
different runs described in the text and in the caption for Table 1-4. The integrated (0-1 00m) biomass is 
also tabulated for each run. The last column shows the corresponding observed phytoplankton nitrogen 
production at BATS, and phytoplankton biomass derived from the observed chlorophyll concentration. 
The model phytoplankton balance does not include nitrogen-fixing phytoplankton (equation 20). All units 
are in mol N m' 2 yr" 1 , except for the integrated biomass, which is in mol N m" 2 . 



Run 1 

Run 2 

Run 3 

Run 4 

Run 5 

Data 

Phyto N Prod 


1.844 

0.873 

1.381 


1.874* 

Grazing 

-0.551 

-1.192 

-0.331 

-0.659 

-0.250 




-0.271 

-0.246 

-0.274 

-0.224 


Mortality 


-0.271 

-0.246 


-0.224 


Advection 





KHHHI 


Mixing 

eehb 


bebi 




Biomass 


14.8 

13.5 


lEBHHHi 

13.6** 


* Phytoplankton nitrogen production is calculated from observed primary production (PP in mg C m" 3 d 1 ): 

0 365 ^ 

Phyto N Prod = — — — J PP dz , where C:N = 6.625 is the Redfield ratio and 12 is the molecular weight 
12 C.N _ l00 

of carbon. The coefficient 0.365 converts mmol N m" 2 d" 1 to mol N m" 2 yr" 1 . 

** Biomass is calculated from observed chl-a (mg m" 3 ) and and chl:N ratio (mg chl-a/mmol N) from model 
(Run 2): 


p = I 


chl-a 

chl:N 


dz 
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Model Validation and Comparison with Previous Model 

With a few exceptions, the model results (run 2) match the observations very well. 
Figure I- 10 shows a comparison between the predicted (solid line) seasonal profiles of 
temperature and salinity with corresponding observed values (black dots). The seasonal 
thermocline, the surface temperature, and surface salinity, are accurately represented by 
the model. Seasonal profiles of the predicted biogeochemical parameters are shown in 
Figures 1-7 and 1-1 1. The overall agreement with the data is good with the exception of 
the summer-fall primary production, which is underestimated by the model, and, as 
explained in the previous section, a much sharper nutricline when compared to the 
observed nitrate and phosphate profiles. The seasonal profiles of TCO 2 , O 2 , DON and 
DOC compare well with observations, except for the winter DON, which is 
underestimated by the model. 

Time series of model (red lines) and observed (green lines for in situ data and blue lines 
for satellite data) chlorophyll concentration and primary production are shown in Figure 
1-12. The model series is daily, while the in situ data were obtained on a nearly monthly 
basis. The satellite data are from monthly composites. The timing of the spring bloom is 
well represented by the model. Again, it is quiet clear that the model under-predicts the 
summer-fall primary production. The predicted chlorophyll concentration tracks the 
observed values well, except for the 1995 spring transition where the fluorometric 
chlorophyll measurements have a few high values (as high as nearly 0.7mg m' 3 ). The 
satellite-derived (OCTS and SeaWiFS) chlorophyll concentrations are from monthly 
composites and therefore are much smoother but still in good agreement with the in situ 
data. Note that the OCTS and model records in early 1997 show anomalous high values 
compared to in situ data. Otherwise, monthly mean satellite values during peaks are low. 

Figure 1-13 shows time series of model (solid) and predicted (dashed) SST and surface 
TCO 2 , calculated alkalinity (regressed on salinity), and predicted ApC02 and sea-air CO 2 
flux. The surface alkalinity concentrations calculated from the model salinity (solid) and 
from the observed salinity (dashed) are very close because the model surface salinity is 
strongly relaxed to observed values. The model (solid) and observed (dashed) TCO 2 
trends are very similar (1.3 and 1.2 pmol kg' 1 yr' 1 , respectively). The predicted seasonal 
and interannual variability in TCO 2 concentration compares well with the observed 
values, except when the predicted SST differs from observations. This is evident during 
the 1995, 1996, and 1997 winters, during which the predicted SST is warmer than 
observations. Also, notice that the model TCO 2 is biased low when compared to the data 
(-3 pmol kg' 1 ). This bias is a result of a corresponding bias in the model SST (+0.48°C), 
probably due to a bias in the total surface heat flux. The time series of predicted sea-air 
CO 2 flux shows larger ingassing during winter than outgassing during summer, despite a 
nearly symmetric summer-winter ApCCh. This can be explained by the larger wind 
speeds attained during winter months when the sea-air ApCC >2 reaches its maximum 
negative value and the sea-air flux reaches its peak ingassing. 


30 



Depth (m) Depth (m) 



Figure I- 10. Seasonal profiles of model (solid lines) and observed (black dots) temperature and 
salinity. 
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Figure 1-1 1. Seasonal profiles of TCO 2 , O 2 , DON, and DOC. The solid lines are model 
profiles (run 2), while the black dots and dashed lines are observed profiles. 
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Figure 1-13. Time series of SST, alkalinity (TA), TCO 2 , ApCC^, and sea-air CO 2 flux (FC0 2 ). The model 
time series are shown as solid lines whereas the observed time series are shown as dashed lines. The TA 
dashed line represents a regression to salinity, not the observed values which are not available. 
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Figure 1-14 shows the predicted (solid line) and measured (black dots) ApC02 for years 
1994 and 1995. Predictions from both the previous (bottom) and current (top) models are 
shown. Some degree of improvement seems to be achieved with the current model, 
especially during 1995. Table 1-6 summarizes seasonal and climatological comparisons 
between predicted and observed biogeochemical parameters. Predicted values are 
tabulated for the previous and current models. The values labeled with asterisks denote 
improved matches between observations and the current model predictions. The 
improvements vary with parameter and season. There were improvements in the model 
simulation of chlorophyll, PAR, nitrate, DOC, primary production, and oxygen, but most 
importantly, the current model achieved good agreement with the data with much more 
realistic background eddy diffusivity. The current model used a constant value of 5x1 O' 5 
mV, while the previous version used variable background eddy diffusivity ranging from 
5xlO' 5 mV to 58xlO* 5 mV. 


Table 1-6. Summary of seasonal and climatological (1992-1998) comparisons between predicted and 
observed biogeochemical parameters. The second line for each parameter tabulates the results from the 
new model. The asterisks highlight values from the new model for which noticeable improvement was 
achieved. 



Winter 

Spring 

Summer 

Fall 

Year 


Model 

Data 

Model 

Data 

Model 

Data 

Model 

Data 

Model 

Data 

PAR 

46.0 

49.3 

83.5 

89.8 

96.3 

99.2 

61.5 

68.9 

71.8 

76.8 


46.4 


86.5* 


99.8* 


63.1 


74.0* 


Chl-a 

0.14 

0.16 

0.14 

0.11 

0.11 

0.04 

0.12 

0.07 

0.13 

0.10 


0.16* 


0.14 


0.11 


0.12 


0.13 


N0 3 

0.43 

0.26 

0.64 

0.16 

0.16 

0.06 

0.09 

0.16 

0.33 

0.16 


0.51 


0.33* 


0.06* 


0T5* 


0.26* 


P0 4 

0.012 

0.009 

0.017 

0.003 

0.002 

0.010 

0.003 

0.008 

0.009 

0.008 


0.002 


0.002* 


0.000 


0.000 


0.001 


PP 

456 

432 

657 

593 

440 

382 

375 

387 

482 

449 


550 


617* 


411* 


218 


449* 


DON 

4.5 

4.4 

4.4 

4.4 

4.5 

4.4 

4.5 

4.3 

4.5 

4.4 


3.7 


4.2 


4.2 


4.3* 


4.1 


DOC 

61.5 

60.9 

60.6 

61.1 

61.8 

63.3 

62.2 

63.5 

61.5 

62.2 


58.0 


57.2 


60.0 


61.5 


59.2 


tco 2 

2049 

2046 

2058 

2055 

2037 

2036 

2030 

2024 

2044 

2040 


2053 


2056* 


2038 


2039 


2047 


o 2 

223 

218 

224 

224 

207 

206 

208 

207 

216 

214 


220* 


224 


208 


207* 


215* 



Nitrate, Phosphate, DON, DOP (mmol m' 3 ), and DOC (yimol kg ') were averaged from the surface to 50m. 
The Chl-n (mg m' J ), TC0 2 (urnol kg' 1 ), 0 2 (pmol kg' 1 ) concentrations are surface values. The observed 
surface irradiance (PAR, in W m' 2 ) is obtained from 1992-1997 averages. The primary production (PP in 
mg C m' 2 d' 1 ) is integrated from surface to 150m for the period of 1992-1997. 


Carbon Budget 

Figure 1-15 shows the mean carbon balance for the BATS region based on the 7-year 
(1992-1998) model simulation. The two boxes in Figure 1-15 represent the dissolved 
inorganic carbon (DIC) in surface water (solubility pump), and the biological pump, 
which takes up CO 2 from the DIC pool via photosynthesis and exports POC and DOC as 
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Figure 1-14. Time series of model (solid lines) and observed (black dots) ApC0 2 . The top tier show 
the series from the current model, whereas the bottom tier show the series from the previous model. 
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Figure 1-15: Averaged (1992-1998) carbon balance for the BATS study area. Fluxes are in mol C m" 2 yr' 
The dilution term is effectively a result of freshwater flux (precipitation minus evaporation). 









9 1 

byproducts. The DOC is remineralized into DIC at a rate of 0.4 mol C m' yr' . The 
biological DIC consumption, i.e., the net community production (2.0 mol C m' 2 yr' 1 ), is 
balanced by the air-sea CO 2 flux (0.3 mol C m’ 2 yr' 1 ), DIC dilution (-0.1 mol C m' 2 yr' 1 ), 
remineralization of DOC (0.4 mol C m' 2 yr' 1 ), and the DIC flux (1.4 mol C m' 2 yr' 1 ) at the 
bottom boundary (350 m). In terms of overall carbon balance, the sedimentation rate 
(POC export), 1 .6 mol C m' yr’ , is balanced by the air-sea CO 2 flux and the DIC 
physical transport through the bottom boundary. 

Conclusions 

The use of the TCMLM contributed to some improvements in the model simulation of 
chlorophyll, PAR, nitrate, phosphate, and oxygen, but most importantly, the current 
model achieved good agreement with the data with much more realistic background eddy 
diffusivity. However, sensitivity tests showed that the TCMLM does not provide 
adequate vertical mixing (convective mixing is ~50m shallower than indicated by the 
data) in winter-spring. In addition, off-line calculations of horizontal divergence of 
nutrients, heat, salt, and dissolved inorganic carbon, using seasonal climatologies, 
revealed that the BATS site lies in a region where horizontal transport of physical and 
biogeochemical properties by the ambient flow is an important process that cannot be 
neglected in studies of nutrient, biomass, and carbon balances. Also, mesoscale eddies 
play an important role in the vertical transport of nutrients (eddy pumping) to the 
euphotic zone, as shown by several previous studies in the region [McGillicuddy et al., 
1999; McGillicuddy et al., 1998b; Siegel et al., 1999]. The role of these eddies on the 
horizontal transport of nutrients, and other biogeochemical properties, remains 
unquantified. Therefore, future studies in the BATS region will require comprehensive 
three-dimensional field studies, combined with three-dimensional eddy resolving 
numerical experiments, to adequately quantify the impact of the local and remote forcing 
on ecosystem dynamics and carbon cycling. 
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Chapter II 


Sensitivity of Global Sea-Air CO2 Flux to Gas Transfer Algorithms, 
Climatological Wind Speeds, Sea Surface Temperature, and Salinity 

11.1. Background 

The global sea-air flux of CO 2 is strongly dependent on the sea-air pCOi difference 
(ApC02) and the gas exchange coefficient. These two components of the air-sea flux 
have similar ranges of variability in the ocean and, therefore, are equally important in 
determining the accuracy of the estimated flux. The gas exchange coefficient is primarily 
a function of the wind speed [Liss and Merlivat, 1986; Tans et al ., 1990; Wanninkhof, 
1992], but the ocean pCOi is highly dependent on temperature and salinity. The ocean 
pCC >2 is a function of the total carbon dioxide (TCO 2 ) concentration, the total alkalinity 
(TA), and dissociation constants for the chemical components of the CO 2 system. 
[Dickson and Miller o, 1987; Signorini and McClain, 2002; Signor ini et al., 2001a; 
Signor ini et al., 2001b]. The dissociation constants are functions of temperature and 
salinity [Dickson and Millero, 1987]. The distribution of surface alkalinity in the open 
ocean is mainly controlled by the factors that govern salinity [Millero et al., 1998]. 
Therefore, accurate measurements of seasonal and interannual changes of surface salinity 
and temperature of the world's oceans are central to studies of global carbon flux. The 
biological uptake of CO 2 is also an important component of the carbon cycle in the 
oceans. However, in this study we do not explicitly address the effect of biological 
uptake of CO 2 . It is assumed that the biological effect is implicitly included in the 
climatological TCO 2 seasonal fields used here. 

11.2. Methodology and Data Sources 

This study is a revision of the approach used by Signorini and McClain [2002]. There 
were two major changes from the previous study. The Wanninkhof [ 1992] gas transfer 
algorithm now uses the 0.39 coefficient for climatological winds instead of the 0.31 
coefficient for short winds. Now the Wanninkhof [1992] and Tans et al. [1990] 
algorithms yield very similar results. Also, the DIC values are no longer kept constant 
when changes of SST and SSS are made. The algorithm of Lee et al. [2000] is used for 
calculating changes in DIC due to variations in SST and SSS. This modification on the 
computational strategy resulted in more internally consistent analyses. 

We used an ocean /?CC >2 model [Signorini et al., 2001a; Signorini et al., 2001b], 
combined with algorithms for TA, TCO 2 , and gas transfer, to estimate the seasonal sea- 
air CO 2 flux for the global oceans. The required input data for the ocean pCC >2 model are 
SST, SSS, TA, and TCO 2 . The normalized seasonal alkalinity (NTA) was obtained from 
the algorithm of Millero et al. [1998], which is a quadratic function of SST. The salinity 
effect on alkalinity is incorporated via dilution (TA=SSSxNTA/35). The regression 
coefficients and reference temperature vary according to specific ocean basins and 
regional domains [Millero et al., 1998]. The seasonal SST and SSS global fields are from 
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the World Atlas 1998 [Conkright el al., 1998]. For the purposes of this analysis the TA 
estimates are assumed to be error-free. 

Existing TCO 2 algorithms are still not accurate enough to retrieve required pC O 2 values 
within reasonable limits. Over the global oceans, ApCOi vary from -100 patm to +100 
patm, (Figure II- 1) but its global mean is always between 0 and -6 patm. An error of 1 
pmol/kg in TCO 2 translates into an error of ~1.5 patm on pQOi. Lee et al. [2000] 
compared TCO 2 fields obtained with their SST and SSS algorithm to the TCO 2 fields 
obtained from the combination of Takahashi et al. (1997) pCC >2 data and TA fields. The 
differences are very large (±40 pmol/kg). For this reason we computed the TCO 2 
seasonal fields from available ocean pC O 2 data [Takahashi et al., 1997], seasonal TA 
fields derived from the Millero et al. [1998] algorithm, and SST and SSS from World 
Atlas 1998 [ Conkright et al., 1998] climatology. Therefore, the TCO 2 seasonal fields 
calculated this way, combined with the TA, SST, and SSS seasonal fields, reproduce the 
seasonal pC O 2 data from Takahashi et al. [1997] exactly. The combination of these input 
data sets was used to perform the sensitivity analysis in this study. The TCO 2 sensitivity 
to SST and SSS is obtained using the Lee et al. [2000] algorithm. 

The flux of CO 2 between the atmosphere and ocean is commonly estimated using a bulk 
formula that relates the flux to the air-sea concentration difference via a gas transfer 
velocity [Wanninkhof, 1992]: 


F co , = ksApC0 2 


(52) 


where the solubility of CO 2 in seawater ( 5 ) is included because the difference (ApC0 2 ) 
is expressed in terms of partial pressures. 

Air-sea gas transfer is typically parameterized in terms of a transfer velocity and the gas 
concentration difference across the air-sea interface. Transfer velocities (k66o) are 
commonly normalized to a nominal Schmidt number, Sc=660, through the following 
relationship which accounts for differences in gas diffusivity 


^660 - k 


660 ' 
Sc . 


- 1/2 


(53) 


Figure II-2 shows the normalized transfer velocity (cm hr' 1 ) as a function of wind speed 
(m s' 1 ) for six different gas transfer algorithms [Liss and Merlivat, 1986; McGillis et al., 
2001; Nightingale et al., 2000; Wanninkhof, 1992; Wanninkhof and McGillis, 1999]. The 
algorithm of Tans et al. [1990], used in this study, is not represented in Figure II-2 
because it is not explicitly expressed as a function of transfer velocity; it is expressed as a 
constant times the wind speed. These six algorithms have quite different behaviors, 
especially above wind speeds of 8 m s' 1 where they diverge significantly. An in situ 
evaluation of air-sea gas exchange parameterizations, using conservative and volatile 
tracers in the southern North Sea, is provided by Nightingale et al. [2000]. They 
concluded that the optimal gas exchange relationship (see Figure II-2) shows a 
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Figure II- 1 . Global map of annual sea-air ApCC >2 from Takahashi et al. (1997) showing the six sub- 
regions chosen for analysis. 
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dependence on wind speed intermediate between those of Liss and Merlivat [1986] and 
Wanninkhof [ 1 992], The Wanninkhof [ 1 992] algorithm is in agreement with the global 
mean transfer velocity determined from the oceanic uptake of bomb-derived radiocarbon 
( 14 C) [Broecker et al., 1985]. Similarly, Tans et al. [1990] adjusted the radon-derived 
parameterization of [Smethie et al., 1985] to fit through the 4 C value. The Wanninkhof 
[1992] algorithm is also in agreement with the gas exchange study of Keeling et al. 
[1998], which is based on seasonal variations in the atmospheric oxygen/nitrogen ratio. 
The gas exchange study of Gabric et al. [1995], based on mass budgeting of 
dimethylsulfide, favors the Liss and Merlivat [1986] algorithm. More recently, McGillis 
et al. [2001] conducted direct co- variance air-sea CO 2 flux measurements in the North 
Atlantic within a significant CO 2 sink. These direct estimates are also in agreement with 
the traditional geochemical isotope constraints. However, McGillis et al. [2001] show 
that measurements for winds in excess of 1 1 m s' 1 indicate a general enhancement of gas 
transfer velocity over previous indirect measurements. They believe that this 
enhancement can be explained by the fact that indirect methods cannot discriminate 
surface process variability such as atmospheric stability, upper-ocean mixing, wave age, 
wave breaking, or surface films. 

The air-sea CO 2 flux (FCO 2 ) is calculated in this study using the following four gas 
transfer algorithms: (I) Liss and Merlivat [1986], Wanninkhof [ 1992], Tans et al. [1990], 
and Wanninkhof and McGillis [1999]. These four algorithms were chosen because their 
behaviors span a wide range of gas transfer velocity values and therefore are suitable for 
the wind speed sensitivity analysis. The seasonal global distribution of atmospheric 
pC0 2 for 1990 was obtained from concurrently available data [Conway and Tans, 1996] 
at three sampling sites, Cold Bay (55°N), Mauna Loa (19.5°N), and the South Pole 
(90°S). The year 1 990 was chosen because it is the reference year used by Takahashi et 
al. [1997] to correct the /CO 2 observations used in their global carbon flux estimates. 
Therefore, our carbon flux estimates are also for 1990. A latitude-dependent global 
seasonal distribution of atmospheric PCO 2 at one-degree resolution was constructed by 
fitting a spline to the three concurrently available records from the three chosen sites. 

Each of the four algorithms were used with three different seasonal wind speed 
climatologies: Esbensen and Kushnir [1981], wind scatterometer data from European 
Remote Sensing (ERS) Satellites, and winds derived from the Special Sensor Microwave/ 
Imager (SSM/I). The Esbensen and Kushnir [1981] climatology was based on data from 
ships of opportunity and was selected for comparison with previous studies that used it. 
Sensitivity tests were also conducted with different SST and SSS climatologies. These 
included the World Ocean Atlas 1998 (WOA98) [Conkright et al., 1998], the World 
Ocean Atlas 1994 (WOA94) [Levitus and Boyer, 1994], and the Reynolds and Smith SST 
[Reynolds and Smith, 1994] data. In order to simulate random error of SST and SSS 
fields, we mixed the WOA94 and WOA98 SST and SSS climatologies in our analyses. 

II.3 Sensitivity to Gas Transfer Algorithms and Wind Speed Climatologies 

Results from the application of the three wind speed climatologies and four different gas 
transfer algorithms are provided for the global oceans and regional domains (Figure II-3) 
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in Table 1. Figure II-3 shows the spatial variability of the annual sea-air CO 2 flux for the 
four gas transfer algorithms and the SSM/I wind speed climatology. Figure II-4 shows 
the seasonal variability of sea-air CO 2 flux for the global oceans and for 6 regional 
oceanic domains derived from the four gas transfer algorithms and SSM/I winds. A 
summary of globally and regionally averaged wind speeds for the three wind data sets is 
provided in Table 2. The sea-air flux varies with the gas transfer algorithm and wind 
speed climatology used. The smallest global sea-air CO 2 flux (Table 1) is obtained with 
the Liss and Merlivat [1986] gas transfer algorithm for all three wind speed 
climatologies. The SSM/I winds produce the lowest global flux (-0.57 Gt/yr), while the 
ERS winds produce the highest flux (-0.70 Gt/yr). The winds from Esbensen and 
Kushnir [1981] yield a global sea-air flux (-0.65 Gt/yr) between those derived from the 
ERS and SSM/I winds. A negative sea-air flux represents uptake of CO 2 by the ocean. 
The largest global sea-air CO 2 flux is obtained with the Wanninkhof and McGillis [1999] 
gas transfer algorithm. This algorithm yields yearly global sea-air CO 2 fluxes of -1.73, 
-2.01, and -2.27 Gt/yr using the SSM/I, Esbensen and Kushnir [1981], and ERS wind 
climatologies, respectively. The Southern Ocean (-0.58 to -1.90 Gt/yr), North Atlantic 
(-0.27 to -0.73 Gt/yr), and North Pacific (-0. 17 to -0.48 Gt/yr) are the largest sinks of CO 2 
throughout the year. As shown in Figure II-4, the North Pacific becomes a weak source 
(0.02 to 0.13 Gt/yr) of CO 2 during the Boreal summer (July- August). The equatorial 
Pacific is the largest source of CO 2 (0.28 to 0.81 Gt/yr) with positive values throughout 
the year. The largest flux in the equatorial Pacific is obtained with the Tans et al. [1990] 
algorithm and SSM/I winds. The Indian Ocean and the equatorial Atlantic are weak 
sources of CO 2 throughout the year (0.04 to 0.16 Gt/yr). 

As shown in Figure II-4, the Southern Ocean has the largest seasonal cycle with a peak in 
ocean uptake during January-February (-1.1 To 2.9 Gt/yr). The North Atlantic and North 
Pacific Oceans have a reduced seasonal cycle, with strongest uptake during fall-winter 
and nearly neutral conditions during the summer months. These regions have small 
outgassing (North Pacific) or small ingassing (North Atlantic) of CO 2 during the summer 
months because the potentially high pC02 values due to elevated SST conditions are 
counterbalanced by an increased pC O 2 drawdown by biological activity. The Indian 
Ocean, equatorial Pacific, and equatorial Atlantic have no distinctive seasonal cycle. 

Our CO 2 flux estimates agree reasonably well with the estimates of Takahashi et al. 
[1997] using Esbensen and Kushnir [1981] winds. Their estimates using the Liss and 
Merlivat [1986], Wanninkhof [ 1 992] , and Tans et al. [1990] gas transfer algorithms are 
-0.71, -1.14, and -1.34 Gt/yr, respectively. Our equivalent estimates are -0.65, -1.28, and 
-1.22 Gt/yr. Lefevre et al. [1999] compiled quarterly maps of A/ 7 CO 2 interpolated from 
pC0 2 measurements in the North Atlantic and the North Pacific Oceans. Their estimates 
of the CO 2 flux for the Northern Hemisphere north of 10°N using SSM/I winds and the 
Liss and Merlivat [1986], Wanninkhof [1992], and Tans et al. [1990] algorithms are 
-0.45, -0.70, and -0.86 Gt/yr, respectively. Our equivalent estimates are -0.52, -1.02, and 
-1.08 Gt/yr. Lefevre et al. [1999] estimates using Esbensen and Kushnir [1981] winds 
and the same three gas transfer algorithms are -0.37, -0.60, and -0.76 Gt/yr, respectively. 
Our equivalent estimates are -0.50, -0.98, -1.03 Gt/yr. Overall, our estimates are 15% to 
35% higher than those of Lefevre et al. [1999]. Since we used pC02 data based on 
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Figure II-3. Annual sea-air CO 2 flux for the global oceans derived from four different gas 
transfer algorithms and SSM/I seasonal wind climatologies. 








Figure II-4. Seasonal variability of sea-air C02 flux derived with SSM/I winds for the global 
oceans and xix different oceanic regions. The four separate lines, solid, dashed, dotted, and 
thick solid, represent the Liss and Merlivat (1986), Wanninkhof (1992), Tans et al. (1990), and 
Wanninkhof and McGillis (1999) algorithms, respectively. 
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Takahashi et al. [1997], and Lejevre et al. [1999] used revised pCOz maps based on 
additional observations, we attribute our CO 2 flux overestimates to the different pCC >2 
data sets used for the estimates. 

11.4 Sensitivity to SST and SSS Climatologies 

Variability in SST and SSS global distributions affect the calculation of TA [Millero et 
al., 1998], TCO 2 [Lee et al., 2000], CO 2 solubility, the Schmidt number, and the ocean 
pC O 2 via the dissociation constants for carbonic acid, boric acid, and sea water, which 
combine to affect the sea-air CO 2 flux. The WOA98 SST and SSS data are used as a 
reference calculation to be compared with the results of other sea-air flux calculations 
using different combinations of SST and SSS climatologies. The Tans et al. [1990] gas 
transfer algorithm was used in all cases. The results are summarized in Table 3, which 
tabulates the FCO 2 means and seasonal ranges, and the mean differences and standard 
deviations of SST, SSS, TA, and TCO 2 . Globally, the mean ocean uptake of CO 2 is 
15.6% larger when the Reynolds and Smith [1994] SST data are used, and the seasonal 
range is 5.5% larger. The SST, TA, and TCO 2 differences associated with the global sea- 
air flux change are -0.1610.53 °C, +0.2611.61 pmol kg' 1 , and +0.86+3.95 pmol kg' 1 , 
respectively. The use of both SST and SSS from the WOA94 data causes an increase in 
ocean uptake of 4.6%, and a seasonal range increase of 5.1%. The combination of 
WOA98 SST and WOA94 SSS causes an increase of 6.4% in CO 2 uptake, with a slight 
increase of seasonal range (1.4%). The global CO 2 flux is least affected (-1.8%) when 
the combination of WOA98 SSS and WOA94 SST is used. The mean deviations in TA 
and TC02 are also small in this case (-0.1111.53 pmol kg' 1 and -0.1013.56 pmol kg' 1 , 
respectively). 

11.5 Conclusions 

Sensitivity analyses of sea-air CO 2 flux to wind speed climatologies, gas transfer 
algorithms, SSS and SST were conducted for the global oceans and regional domains. 
Large uncertainties in the global sea-air flux are identified, primarily due to the different 
gas transfer algorithms used. However, two of the algorithms provide very close results. 
Globally, there is less than 9% difference between the Tans et al. [1990] and the 
Wanninkhof[ 1992] algorithms, which produce sea-air fluxes within the mid-range of our 
estimates. The use of three different wind climatologies with these two algorithms 
produces uncertainties of 7 to 10%. The other two algorithms [LAs and Merlivat, 1986; 
Wanninkhof and McGillis, 1999] differ by a factor of two and produce sea-air fluxes at 
the low end [Liss and Merlivat, 1986] and high end [ Wanninkhof and McGillis, 1999] of 
our estimates. The sensitivity of the sea-air flux to SST and SSS is similar in magnitude 
to the effect of using different wind climatologies. Globally, the mean ocean uptake of 
CO 2 changes by 5 to 16%, depending upon the combination of SST and SSS used. 
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